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ABSTRACT 

Ionizing stellar photons heat the upper regions of planetary atmospheres, driving atmospheric mass 
loss. Gas escaping from several hot, hydrogen-rich planets has been detected using UV and X- 
ray transmission spectroscopy. Because these planets are tidally locked, and thus asymmetrically 
irradiated, escaping gas is unlikely to be spherically symmetric. In this paper, we focus on the 
effects of asymmetric heating on local outflow structure. We use the Athena code for hydrodynamics 
to produce 3D simulations of hot Jupiter mass loss that jointly model wind launching and stellar 
heating via photoionization. Our fiducial planet is an inflated, hot Jupiter with radius = 2.14i?jup 
and mass Mp = 0.53Mjup. We irradiate the initially neutral, atomic hydrogen atmosphere with 
13.6 eV photons and compute the outflow’s ionization structure. There are clear asymmetries in 
the atmospheric outflow, including a neutral shadow on the planet’s nightside. Given an incident 
ionizing UV flux comparable to that of the Sun, we find a steady-state mass loss rate of ^2 x 10^^ 
g s“^. The total mass loss rate and the outflow substructure along the substellar ray show good 
agreement with earlier ID models, for two different fluxes. Our 3D data cube can be used to generate 
the outflow’s extinction spectrum during transit. As a proof of concept, we find absorption of stellar 
Lya at Doppler-shifted velocities of up to ±50 km s“^. Our work provides a starting point for further 
3D models that can be used to predict observable signatures of hot Jupiter mass loss. 

Subject headings: hydrodynamics - planets and satellites: atmospheres, gaseous planets - planet-star 
interactions 


1. INTRODUCTION 

Photoionization by high energy stellar radiation heats 
the upper layers of planetary atmospheres. This heating 
drives thermal mass loss and can thus play a substantial 
role in a planet’s atmospheric evolution. Strongly irradi¬ 
ated hydrogen-rich planets are most susceptible to mass 
loss, and transit observations of close-in giant planets at 
UV and X-ray wavelengths have revealed atmospheric 
escape. However, models are required to translate these 
observations into constraints on mass loss rates or out¬ 
flow structure. While models have been developed for hot 
Jupiter mass loss, none have yet consistently modeled the 
heating and three-dimensional (3D) structure of atmo¬ 
spheric escape. The inherent asymmetry in the physics 
of atmospheric escape, especially asymmetric irradiation 
expected of tidally locked hot Jupiters, necessitates full 
3D modeling. To examine how asymmetric heating af¬ 
fects the structure of the outflow near the planet, we 
develop a 3D, self-consistent model of mass loss driven 
by photoionization heating. 

The first indication of hot Jupiter mass loss came 
fro m the transmission spectroscopy of 
al. (2003), who observed decreased stel 


Vidal-Madjar et 
ar Lya emission 
Hurmg transits of the hot Jupiter, HD 209458b. The 
depth of the transit in Lya, 15% ± 4%, was ten times 
larger than the optical transit depth of 1.5% (Gharbon- 
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neau et al. 2000 Henry et al.|2000), suggesting absorp- 
tion by a neutral hydrogen atmosphere larger than the 
planet’s Roche lobe. Absorption out to Doppler equiva¬ 
lent velocities of ±100 km s“^ indicated that the neutral 
gas either moved at high velocities or had a large column 
depth, enhancing the line’s naturally broadened wings. 
Subsequent studies have confirmed stellar L ya absorp¬ 


tion by ext ended gas around HD 209458b (Eh renreich 
et al. 2QQ8| and the hot J upiter HD 189733b (ju 


^ecave- 


lier des Etangs et al.||201Q), with temporal variations m 


the measured absorp tion ( Ben-Jaffel|2008 Lecavelier des 


Etangs et al. 20 12|). Absorption in CIU Ol, and Silll 
for HD 209458b fVidal-lVladjar et al.||2004l iLinsky et al. 
IT 


but see 


Ballester V Ben-Jaf[:el|2015t and in X-rays 


2010 

for HD 189733b ( [Foppenhaeger et al.|2013[ ) suggest that 
the atmospheric outflows of these planets are metal en¬ 
riched, although interpretation of th e X-ray signal may 


be co mplicated by stellar variability (Llama & Shkolnik 


2015). Metal line absorpti on has also been de tected for 


the hot Jupiter WASP-12b (Haswell et al. 2012[), which i s 
thought to be overflowing its Roche lobe ( ET et al.|2010 ). 
Stellar Lya absorption observations of a transiting hot 
Neptune, GJ 436b, suggest it has an asymmetric out flow 
structure (|Kulow et al.||2014t |Ehrenreich et al.||2015|) . 

Stellar heating induces planetary mass loss through hy¬ 
drodynamic escape. This heating is primarily due to 
extreme UV stellar radiation, which photoionizes neu¬ 
trals in the planet’s upper atmosphere, liberating elec¬ 
trons that heat the gas through collisions. This pho¬ 
toionization heating creates pressure gradients that ac¬ 
celerate the gas from subsonic to supersonic velocities. 
Gonsequently, heated gas moves outward in the form of 
a pla netary wind, similar to the structure of the solar 
wind (Parker| 1958). Gas exceeding the escape velocity 
or displaced beyond the Roche lobe becomes unbound 
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and escapes into space. For the large ionizing fluxes re¬ 
ceived by hot Jupiters, hydrodynamic escape is the most 
efficient mass loss proce ss. Other mass loss mechanisms, 
including Jeans escape ( Chamberlain||1963 ), are less effi¬ 
cient because they operate on individual particles, rather 
than a collective fluid. 

One-dimensional (ID) hydrodynamic models of 
energy-li mited escape from ho t Jup iters were first calcu- 


lated by ham mer et al. 
to explain the | V idal- 


and 


adjar et al. 
Based on work 


Baraffe et al. (2004) 
observations 


done for the early 


of HD 209458t 

Earth and Venus ([Gross 1972 Watson et al. 1981|), 
these pioneering calculations assumed that all energy 
deposited by stellar radiation goes into heating. The 
heat is conducted to lower radii and drives gas to larger 
radii, where it can absorb more flux, thus enhancing 
the outflow. These models predicted catastrophic evap¬ 
oration fo r hot Jupiters. Other energy-limited models 
(including |Hubbard et al. 2007) and more detailed ID 
models which have accounted for chemistry, h eating 
and cooling, tidal gravity, a nd the stellar wind (|Yelle 
2004[ IGarcfa Muhoz| |20Q7| Murray-Clay et ^ |2009[ 
hereafter MG09) have suggested less dramatic mass loss 
for hot Jupiters. Hot Jupiters, unlike the early Earth 
and Venus, do not efficiently conduct heat downward 
(jGarcfa Muhoz||2007| MC09). 

hJot modeled by ID studies are the inherent asymme¬ 
tries in atmospheric escape processes. The stellar wind’s 
pressure confinement, rotation from the Goriolis force, 
and magnetic fields add to the asymmetry of escaping 
atmospheric gas. Day and night differences due to the 
stellar wind have been captured by mass loss models in 


P .. .. 

in 3L) ([Bisikalo et al. 20131 [LJamn^^Vah 120^ 


,__ W _ _ _, 

2D ([Stone fc Proga||2QQ^jTrembhn fc Ghiang||2Q13 ) and 


and other 3D simulations (ISchneiter et a. 


'hese 


Cohen 


et al.|2Qlll IBourrier fc Lecavelier des Etan^ 


sakos et al.|2015 ), which variously include orbital motion, 

radiation pressure, or magnetic fields, have not directly 
included the planetary wind’s production by ionizing ra¬ 
diation from the host star. Instead, they have initial¬ 
ized the temperature at the wind base and used this to 
generate a plane tary wind. The recent work of [Owen] 
& Adams (2014) moves toward a more self-consistent 
picture by simulating ionization-driven winds from hot 
Jupiters, with stellar winds and magnetic fields, in 2D. 
Still, there is a need for 3D simulations which take into 
account photoionization, and this motivates our work. 

As a first step toward realistic models of hot Jupiter 
mass loss, we present a new 3D hydrodynamic model 
of atmospheric escape with self-consistent heating. We 
focus on how asymmetric heating affects the flow near the 
planet. The model and results are organized in this paper 
as follows. In Section^ we describe the physics included 
in our simulation. In Section[3]we describe our simulation 
setup and our initial conditions. In Section]^ we present 
our results of the time-evolved wind structure, mass-loss 
rates, and comparisons to ID models. In Section we 
provide our estimate of the predicted Lya transmission 
spectrum mid-transit. In Section we conclude with a 
summary and a discussion of future extensions. 

2. MODELING HYDRODYNAMIC ESCAPE 

To model hydrodynamic escape we conduct 3D radi¬ 
ation hydrodynamic simulations. Although the upper 


atmosphere of a hot Jupiter is low density, the mean free 
path remains small compared to the scale height, justify¬ 
ing the fluid approximation. To attain sufficient dynamic 
range to resolve the upper atmosphere, we include only 
the planet in our simulation domain. The star resides 
outside of the computational boundary and exerts its in¬ 
fluence through the gravitational potential and through 
ionizing photons entering a single boundary. 

We use th e publicly available grid-based code Athena, 
version 4.1 (Stone et al. 2008), to solve the equations 
of ideal hydrodynamics. We impleme nt an additional 
modul e for photoionizing radiation from jKrumholz et al. 


(2007), as described in Appendix[A| Our modihed version 
of Athena, with initial condition^les, is freely available 
for download and usej^ In this section, we describe the 
numerical implementation of our model, and our initial 
conditions can be found in Section [3l 


2.1. Hydrodynamics 

We use Athena to solve the following set of hydrody¬ 
namic equations, including gravitational, radiative, and 
chemical evolution source terms: 


|+V.(,v)=0. 

(1) 

d 

— (/9v) + V ■ (pvv) + VP=-/)V$, 

(2) 

ri F 

—+V-((E + P)v) = g-£, 

(3) 

^ + V■(p„v)=7^-I, 

(4) 

where p is the total density, pn is the density of neutral 
gas, V is the velocity, P is the thermal pressure, and 4> is 
the gravitational potential. The total energy density E, 
excluding the chemical potential energy, is 

V • V 

E — € + p , 

(5) 


where e is the internal energy density excluding the chem¬ 
ical potential. Omitting the chemical potential allows us 
to use the usual relationship between pressure and inter¬ 
nal energy for an ideal gas. 


and to adopt an adiabatic equation of state with 7 = 5/3 
as a constant, appropriate for either the atomic or ionized 
gas expected to be found in the upper atmosphere of a 
hot Jupiter. Excluding the chemical energy and treating 
7 as constant in this manner is reasonable because we 
are in a regime where collisional ionization is negligible, 
and thus there is a strong separation of scales between 
the mean particle kinetic and chemical energies. In the 
energy equation. Equation Q and jC are the rates of 
radiative heating and cooling, respectively. In the con¬ 
tinuity equation for the neutral density. Equation IZ 
and X are the rates of recombination and ionization, re¬ 
spectively. 

^ https://github.com/tripathi/Atmospheric-Athena 
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2 . 2 . Gravitational potential 

We employ a static gravitatioual poteutial from both 
the plauet aud the host star, iucludiug the coutributiou 
from the ceutrifugal or tidal term: 

GMp GM* lGM*r2 

r r* 2 a3 ’ ^ ^ 

where Mp is the mass of the plauet, is the mass of the 
star, a is the distauce betweeu the ceuters of the plauet 
aud star, r is the local distauce to the plauet, aud is 
the local distauce to the star. 


Our treatmeut of photoiouizatiou aud recombiuatiou 
assumes the case-B couditiou aud the ou-the-spot ap- 
proximatiou, so that the gas is optically thick to iouiz- 
iug photous iu regious where we cousider recombiuatiou. 
For recombiuatious to the grouud level, we assume that 
emitted iouiziug photous will iouize at the same locatiou, 
effectively cauceliug out that recombiuatiou. Thus, the 
time rate of chauge of the recombiued deusity, iu uuits 
of g cm“^ s“^, is giveu by 

'll = , (H) 


2 . 3 . Ionization balance 

As uoted iu Sectiou [2Tj we track both the ueutral aud 
the iouized gas iu the computatioual domaiu. We model 
the chauges iu ueutral deusity based ou rates of photoiou¬ 
izatiou aud recombiuatiou. The former depeuds ou the 
available stellar iouiziug flux, while the latter depeuds ou 
the uumber deusities of ious aud ueutrals iu the gas. 

We cousider a simulatiou box whose distauce from the 
star is siguificautly larger thau the size of the box, /box- 
We ueglect geometric spreadiug of the stellar radiatiou 
field aud simply treat it as a plauar frout of radiatiou 
euteriug the box. Siuce the uear edge of the box is 4.5 
/box from the star aud the far edge is 5.5 /box, the chauge 
iu flux betweeu the uear aud far edges of the box aud 
the ceuter, where we have uormalized, is 20%. We are 
therefore makiug au error of this order by ueglectiug ge¬ 
ometric spreadiug. A plauar frout of radiatiou euters the 
box with photou flux Fq. Iu this approximatiou, the flux 
at a distauce x iuto the box is giveu by 

F(x)=Foe-"W, (8) 

which depeuds ouly ou the optical depth r, defiued as 

r(x) = / nuCTphdi, (9) 

Jo 

where cTph is the photoiouizatiou cross sectiou aud ur is 
the ueutral uumber deusity. 

Giveu the abuudauce of hydrogeu iu hot Jupiter at¬ 
mospheres, we cousider a pure hydrogeu atmosphere so 
that uh = where the meau gas mass per hydro¬ 

geu uucleus is ph = '^h = 1-67 x 10“^^ g aud mn is the 
mass of a hydrogeu atom. For hydrogeu photoiouizatiou 
at the threshold euergy of 13.6 eV, the cross sectiou is 
o’ph = 6.3 X 10“^^ cm^. Iu reality, the cross sectiou will 
depeud ou the photou frequeucy as roughly how¬ 
ever, siuce the chemical state is mostly coutrolled by 
photous uear the iouizatiou threshold for optically thiu 
gas aud the euergetics are coutrolled by photous that 
spau much less thau a factor of 2 iu frequeucy, we simply 
adopt the threshold cross sectiou for all purposes. While 
Equatiouj^is true for auy frequeucy, we assume that our 
iouiziug source is mouochromatic, for comparisou with 
previous work. We review our choice of mouochromatic 
flux iu Sectiou o The calculated flux is used to deter- 
miue the time rate of chauge of the iouized deusity, iu 
uuits of g cm“^ s“^, as 

^ph = Crphpn7^(x). (10) 

The photoiouizatiou rate is simply Xph /Ph • 


where = 2.59 x 1Q~^^(T/1Q ^ K)~^-^ cm^ 
case-B recombiuatiou coefhcieut ( |Osterbrock|1989 ), nH+ 
is the uumber deusity of protons, aud Uq is the num¬ 
ber density of electrons. We assume that = nH+ = 
(p ~ Pn)/PH- Note that the recombination rate is simply 
7 ^/Ph- 

We omit collisional ionization from our treatment be¬ 
cause it is negligible for our problem, since MC09 found 
that temperatures at the atmosphere’s wind base are 
^ 10^ K. The gas temperature T is determined from the 
total energy density E and momentum density p using 


T = 


7-1 


P 


2 p J k 


( 12 ) 


where the mean gas mass, p = xpi + (1 — x)pAr, depends 
on the ionization fraction, x = 1—pnl P, the mean particle 
mass of ionized species pi, and the mean particle mass 
of neutral gas pw- For our hydrogen gas, we use pi = 
mH/2 and pw = = 1.67 x 10“^^ g. We note that 

though the true value of p is p = mH/(l + x) for atomic 
and ionized hydrogen, the above expression, chosen for 
compatibility with previous code development, provides 
a good approximation. 


2 . 4 . Heating and cooling 

Radiation not only ionizes the gas, but it also con¬ 
tributes to its heating and cooling. Each ionizing photon 
imparts its energy in excess of the ionization threshold 
to the newly liberated electron, which then heats the gas 
through collisions. Each photon contributes an energy, 
er, to heating the gas, so that the photoiouizatiou heat¬ 
ing rate, per unit volume, is 


0ph = erCrphUHE(x). 


(13) 


We use au ISM appr opriate heating rate of ep = 2.4 eV 
( Whalen et al.|[2QM ), chosen to validate our ionizing ra- 
diative transfer code ag ainst HII region ionization fronts, 
discussed in Appendix |A.2[ While we use this heat¬ 
ing rate for both our validation tests and our mass loss 
model, we note that for a hot Jupiter atmosphe re absorb¬ 
ing a quiet so lar Lyman continuum spectrum, |TrammeII 
et al. (2011) calculated er = 2.7 eV, assuming 100% 
euergy Mepositiou efficiency. Efficiency values may be 
smaller (Koskinen et al. 2014), making our choice of heat¬ 
ing rate an upper limit tor rnass loss. 

Hydrogen recombination emits photons, which can es¬ 
cape and thus cool the gas. The radiative re combiuatiou 
cooling rate, from a linear fit to values in Osterbrock 
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Fig. 1.— Our computational grid has five levels of resolution 
(each level is in a different color), to resolve the scale height of 
the planet (gray). The grid spans 1.5 X 10^^cm (lOi^p) in each 
dimension, with a highest resolution of 1/128 Rp. This 2D slice in 
the x — y plane shows how the grid is distributed across processors, 
as denoted by the thick lines. 


(1989), is 


(6.11 X 10 


-10 


Q —1 

cm s 




-0.89 


(14) 

where k is the Boltzmann constant. 

We also include the rate for cooling from neutral atoms 
th at are collision ally excited and emit hya photons, given 
by Black| ( |1981 ) as 


^Lya = (7.5 X 10-1® erg 

___ (15) 

Menager et al.|(2013) suggest that this formula may over- 
estimate the Lya cooling rate by up to an order of mag¬ 
nitude for hot Jupiter atmospheres, due to non-LTE and 
other radiative transfer effects. As we discuss in Section 
4.1.2[ Lyof cooling is negligible almost everywhere for our 


planet’s parameters, so such a difference does not change 
our results, but it may affect hotter and denser planets. 


TABLE 1 

Simulation Parameters 


Parameter 

Value 

Planet 


Mass, Mp [g] 

IqSO 

UV photon absorption radius, Rp [cm] 

1.5 X IQio 

Density at i?p, pp [g cm“^] 

0 

1 

CJi 

Isothermal sound speed at i?p, Cs,p [cm s“^] 

3 X 10^ 

Hydrogen photoionization cross section, (jph [cm^] 

6.3 X IQ-i® 

Mean mass per hydrogen nucleus, pn [g] 

1.67 X 10-24 

Star 


Mass [g] 

2.0 X lO^^ 

Orbital distance, a [cm] 

7.48 X lO^i 


origin of a 3D Cartesian grid. The computational volume 
corresponds to a physical size of ( 10 i?p)^, discretized into 
a base grid of 80^ cells. To resolve the minimum atmo¬ 
sphere scale height when the wind is launched, we use 
five levels of grid refinement around the planet, yield¬ 
ing a finest resolution of (1/128 We show in Ap¬ 

pendix that our resolution is sufficient for our results 
to reach numerical convergence. At this resolution, the 
influence of the Cartesian discretization on the spherical 
atmosphere is minimal. 

We work in the frame corotating with the planet’s or¬ 
bit. We assume that the planet is at a fixed distance from 
the star. In this work, we omit the Coriolis force. We do 
not include the rotation of the planet since hot Jupiters 
are tidally locked. We inject stellar UV flux into one face 
{—x) of our box and calculate the optical depth to pho¬ 
toionization from this boundary. For the fluid variables, 
all faces of our box have outflow boundary conditions. 

3.1. Planetary atmosphere 

To study hydrodynamic escape, we retain a large, well- 
resolved atmosphere to serve as a mass reservoir from 
which a wind can be launched. The atmosphere is in 
hydrostatic equilibrium, with a gas density profile: 

dP _ _GMpp 
dr 

and in agreement with our adiabatic equation of state, 
we assume a polytrope 

P = Kp\ (17) 


2.5. Numerieal Algorithm 

We use Athena’s Roe’s linearized Riemann solver, with 
default second-order spatial reconstruction of the fluid 
variables and the directionally n nsplit corner trans port 
upwind (CTU) integrator in 3D^(| Stone et al. [ |2QQ8 ). To 
avoid the carbuncle instability 


uirk 


1994p , which we 


find when photoionization heated gas advects around the 
planet and converges on the night side, we use Athena’s 
built-in H-correction (Stone et al.||2008), to add dissipa¬ 
tion when strong shocks are aligned with the grid. Our 
fiducial simulation length is roughly four orbital periods, 
by which point the wind has reached a steady state. 


3. INITIAL CONDITIONS 

In this section, we describe our choice of initial con¬ 
ditions for the computational domain, planetary atmo¬ 
sphere, and the host star. As shown in Figure we 
initialize our computational domain with a planet at the 


where K is di constant of proportionality and 7 is the 
adiabatic index. Integrating from a radius Rp (with cor¬ 
responding density Pp) to a radius r yields 


Patm {j ') 


/I 

7 K \r 




l/(7-l) 


(18) 

For an ideal gas, K can be determined from the isother¬ 
mal sound speed and p using 


K = 


1 —'T 2 

p 


(19) 


We normalize these parameters at Rp, the radius where 
ionizing photons are absorbed and from which the wind 
is launched. 

Resolving the entire planet is computationally pro¬ 
hibitive and unnecessary, so we create an artificial inner 
boundary for the atmosphere. We set this inner bound¬ 
ary at Tb = 0.75Rp, ^ 4 scale heights below the wind 
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lauuchiug poiut, which is more thau sufficieut to maiu- 
taiu a reservoir of atmospheric gas. luterior to rb, the 
deusity aud temperature profiles are fixed at their hy¬ 
drostatic values at every timestep, but uou-zero fluxes 
are permitted across the bouudary. To preveut the 
deusity from divergiug at the origiu, as expected from 
Equatiou 18, we set p{r < vq) = Patm(^o) at a radius 
vq = O.bRp. Our model is iuseusitive to the choice of tq, 
so loug as it is smaller thau rb by the uumber of cells 
used for the Riemauu solver’s recoustructiou method. 

To maiutaiu the stability of the simulated atmosphere, 
we match the pressure at the atmosphere’s outer edge 
to that of a statiouary aud uuiform ambieut medium. 
Uulike the pressure, the deusity of the ambieut medium 
is discoutiuuous from the plauet. To eusure that the 
backgrouud gas does uot iuflueuce the developmeut of 
wiuds from the plauet, we require it to be at a low euough 
deusity (or equivaleutly, a high euough temperature) so 
that the iuevitable accretiou of material outo the plauet 
has a miuimal effect ou the atmospheric structure aud 
wiud lauuchiug. Combiuiug these coustraiuts, we use 
the followiug profile to specify the deusity of the gas iu 
the eutire domaiu. 


Patm('^o) 

P(r) = { pB.tm{r) 

,Patm(re) ' 10“^ 


r < ro, 

ro<r < Te, (20) 
r > re, 


where is the edge of the plauet. The pressure profile 
is 


(Kp^tmiroV 
P{r) = I Kp^tmiry 
[Kp^tmireV 


r < ro, 
ro < r < re 
r > re- 


( 21 ) 


Siuce we model a spherical plauet ou a Cartesiau grid, the 
pressure is uot perfectly matched betweeu the ambieut 
gas aud the edge of the plauet’s atmosphere. Neverthe¬ 
less, we hud that the chauges iu our atmosphere over a 
souud crossiug time are substautially smaller thau the 
outflows that develop. 

Our fiducial parameters are summarized iu Table 
We model a low mass aud exteuded hot Jupiter with 
Rp = 1.5 X 10^^ cm = 2.14Rjup aud Mp = 10^^ g = 
0.53Mjup. These values are similar to WASP-17b, the 
most iufi ated exoplauet t o-date, with Rp = 1.97 ± 
Q.Q6R.Tiin (|Beuto et al. | 2Q 14[) aud Mp = 0.477±0.033Mjup 
( Southworth et al.|2U12 r~r^ is importaut to uote that the 
radius that we quote as Rp for our model correspouds to 
the absorptiou radius for iouiziug photous, which is larger 
thau the observed optical plauetary radius. Our choice 
of a low mass, low deusity plauet reduces the computa- 
tioual cost to resolve the upper atmosphere aud lauuch 
au outflow. 

Followiug the parameter study of MC09, which fouud 
that mass loss rates are iuseusitive to wiud-base temper¬ 
atures < 10^ K, we set Cs^p = 3 x 10^ cm s“^, or equiva¬ 
leutly Tp = 1.1 X10^ K. Our defiuitiou of Rp requires that 
over a scale height evaluated at Rp, = Cg/(GMp/Rp), 
the optical depth reaches uuity. This correspouds to a 
uumber deusity, Up = 6 x 10^ cm“^. To eusure that the 
gas is optically thick, we exteud the plauet’s atmosphere 
to a radius Vq = 1.02Rp, where p(re) = Pp/10. We begiu 
with a ueutral hydrogeu atmosphere with the deusity at 


^p: Pp = Pu^p = 10“^^ g cm“^. We uote that though 
we set our iuitial couditious so that iuitially r ^ 1 at Rp, 
the gas is allowed to self-cousisteutly choose where the 
T = 1 surface lies. While the plauet begius completely 
ueutral, the backgrouud gas is fully iouized. Startiug 
with a fully iouized aud optically thiu backgrouud allows 
us to track the evolutiou of the plauet aloue. 


3.2. Stellar radiation 

We iuclude the host star’s gravity iu our static gravita- 
tioual poteutial. We set the mass aud radius of the star 
as = Mq aud R^ = Rq. The stellar radius is used 
to calculate the predicted hya extiuctiou duriug trausit. 
Typical of hot Jupiters, the star is located at a distauce 
of 0.05 AU from the plauet. 

Siuce hot Jupiters are tidally locked to their host stars, 
stellar flux is coustautly received by the same face of the 
plauet. We iuclude the stellar flux as a plaue-parallel 
source of iouiziug radiatiou iucideut ou oue si de of our 
computatioual box. As described iu Sectiou |2.3[ the 
plaue-parallel approximatiou is justified by the relative 
sizes aud orbital separatiou of the plauet aud the star. 
We treat the stellar flux as a mouochromatic source, 
rather thau a full stellar spectrum. This choice allows us 
to make direct comparisous with the MC09 ID model, 
which also uses a mouochromatic source. 

We illumiuate the plauet gradually, slowly iucreasiug 
the flux over two orders of maguitude, usiug the fuuctiou. 


for a giveu choice of Fp. The “ramp up” timescale is 
choseu to be siguiflcautly louger thau the advectiou time 
for gas moviug arouud the plauet. By iucreasiug the 
photou flux slowly, we allow the system to gradually relax 
to equilibrium. 

To examiue chauges iu mass loss as a fuuctiou of 
flux, we ruu our simulatiou with two differeut stellar 
flux values. For our fiducial model, we use a maxi- 

14 ^^-2 „-l 
& , 
-2 ^-1 


mum steady state photou flux of 2.02 x 10"" cm " s 


which correspouds to IO.IRq foi* Fp = 2x 10^^ cm“^ s“ 
Siuce thi s value of Rp is com parable to the solar UV 
flux from Woods et al. (1998) scaled to our orbital dis¬ 
tauce of 0.05 AU, our fiducial model represeuts a rel¬ 
atively youuger, hotter star thau the Suu. To study 
a Suu-like star, we also ruu a secoud model, where 


Fq = 2 X 10^^ cm“^ s“^. Because we use siugle fre- 
queucy photous, there is uot a oue-to-oue correspoudeuce 
betweeu our choice of Fp aud the solar flux; a differeut 
frequeucy rauge would uecessitate a differeut Fp. 


4. RESULTS 

The s truc ture of the plauetary outflow is described iu 
Sectiou 4.1 , with wi uds ou the day aud uight sides de¬ 
scribed iu |4.1.1 aud |4.1.2 , respectively. The mass loss 
rate for our hducial mod el, as well as the lower-flux 
model, is preseuted iu |4.2[ Agreemeut bet weeu our sim¬ 
ulatiou aud ID models is discussed iu 14.31 


4.1. Wind structure 

Withiu au orbital period of illumiuatiou, the plauet 
develops a steady-state, trausouic atmospheric outflow. 
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Density [g cm ^] 


5*io-2i io-19 


Neutral fraction 
10-4 10-3 


Temperature [K] 

105 i.sho® 


Line-of-sight X-velocity [km s fl 



-25 -12.5 

6 12.5 25 



tL 

55 


x/Rp 

Fig. 2. — Atmospheric expansion on the dayside, with advection toward the nightside, as seen in this time evolution of density [g cm“^] 
with velocity vectors overplotted (left), neutral fraction (second column), temperature [K] (third), and x-velocity [km s“^] (right), for a 
slice through the midplane. Top row (t=0s): initial conditions, second row (3 X lO^s): the atmosphere heats on the dayside and advects 
to the nightside, third row (6 x lO^s): the outflows are tidally extended along the axis to the star, fourth row (1.85 X lO^s): the outflows 
continue at larger radii, with lower density, and Anally bottom row (1.39 x 10®s): the steady state outflow. The star illuminates the planet 
from the —x boundary, and the mid-transit observer is at -\-x. 
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Fig. 3. — Steady-state outflow radial velocity [km s“^], shown 
with the adiabatic sonic surface (solid) and the escape surface 
(dashed) in the midplane. For reference, the Roche lobe radius 
is 2.7Rp. 

After reachiug steady state at ^ 3 x lO^s, escape coutiu- 
ues iu steady state for the remaiuder of our simulatiou: 
three orbital periods or four souud-crossiug times for our 
simulatiou box. The trausitiou to steady-state cau be 
seeu iu th e tim e evolutiou of the mass loss rate, discussed 
iu Sectiou [4^ 

The wiud’s evolutiou from iuitializatiou to steady state 
is illustrated iu Figure This aud o ther midplaue visu - 
alizatious were geuerated usiug Visit ( [Childs et al.|2Q12 ). 
Siuce our model is axisymmetric, the outflow structure 
is the same iu both the x — z (showu here) aud the x — y 
plaues. The day side exhibits a stroug, iouized outflow. 
The uightside outflow is slightly weaker, with ueutral 
gas iu the plauet’s shadow. We discuss these day-uight 
differeuces aud the uuderl yiug euerg etic aud iouizatiou 
processes iu Sectious [d.l.lj aud [4.1. 2| 

As showu iu Figure [^ the steady state outflow is 
trausouic. Its radial velocity reaches the local es¬ 
cape speed after exceediug the adiabatic souud speed. 
The flow becomes supersouic at distauces from the 
plauet comparable to the Roche lobe radius, RRoche = 
[Mp/(3M*)]V3a = 4.1 x 10i°cin = 2.7Rp. The escape 
surface is largely outside of the Roche lobe radius aud is 
closest to the plauet aloug the substellar ray. 

Tidal gravity is respousible for elougatiug the outflow 
aud coutributes to the asymmetric escape surface. Cou- 
sequeutly, the velocity field of the outflow observed dur- 
iug a plauetary trausit should be domiuated by the liue- 
of-sight velocity. 

Throughout the results sectiou, we discuss the wiud 
lauuchiug iu time, so as to uuderstaud its flual structure 
aud dyuamics. Because our iuitial couditious do uot re¬ 
flect the properties of a uewly formed hot Jupiter, this 
time evolutiou does uot represeut the early evolutiou of 
a hot Jupiter’s outflow. However, the model’s time evo¬ 
lutiou may be relevaut for plauets orbitiug very active, 
flariug stars, whose flare timescales cau be comparable 
to our simulated flux ramp up time. 

4.1.1. Day side flow: directly launched hy photoionization 

heating 


Temperature 


- 1000 


5000 

10000 

23000 


x/Rp 

Fig. 4.— The temperature [K] of the steady-state outflow ex- 
hibits a day-night asymmetry, seen in this midplane slice. While 
dayside temperatnres remain < 10^ K, nentral nightside gas ap¬ 
proaches 10‘^K. 

Ou the plauet’s dayside, photoiouizatiou heats the at¬ 
mosphere to ^ 7000K, by depositiug euergy at Rp, where 
r = 1 to photoiouizatiou. Irradiated gas above ac¬ 
celerates aud moves outwards. This outward expausiou, 
hereafter referred to as PdV work, is the primary coolaut 
of the gas. Radiative recombiuatiou cooliug is au order of 
maguitude smaller. Because our plauet has low surface 
gravity, its outflow does uot achieve the lO^K temper¬ 
ature required for substautial hya cooliug, as showu iu 
Figure Lower temperatures are sufhcieut to accelerate 
gas to the plauet’s escape speed. 

To assess iu more detail the relative coutributious of 
various heatiug, cooliug, aud iouizatiou processes to the 
wiud’s structure, we recast Equatious aud as 

P 

^ + — (V • Vp) + 0ph — C^ec ~ >^Lya = 0 
^ (23) 

—V • Vx + —— ^ = 0. (24) 

/i ph Ph 

The first two terms iu Equatiou [^ represeut the chauge 
iu iuterual euergy aud the PdV work, respectively. The 
first term iu Equatiou represeuts the advectiou of ious 
out of a giveu cell. 

Above the wiud base, PdV cooliug aud photoiouiza¬ 
tiou heatiug coutribute to euergy balauce, as showu iu 
Eigure which displays the terms of Equatiou aud 
iu steady state. Near the wiud base, heat is stored 
iu iuterual euergy, as showu by the uegative iuterual eu¬ 
ergy below l.lRp aloug the substellar ray. Earther from 
the wiud base, the outflow is driveu by local photoiou¬ 
izatiou heatiug, as well as the heat that was deposited 
by photoiouizatiou lower iu the atmosphere aud stored 
iu iuterual euergy. 

Due to the large photoiouizatiou rate of ^ 
10^ cm“^ s“^, the dayside outflow is everywhere iouized. 
The dayside is iu iouizatiou balauce above the wiud base, 
as seeu iu the lower pauel of Eigure Near the wiud- 
base, photoiouizatiou coutributes a steady source of ious, 
which are advected away. lou advectiou is ouly impor- 
taut uear the wiud base, where the fractiou of the gas 
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Fig. 5. — Steady-state energy (up per) and ionization (lower) balance along the substellar ray on the dayside (left) and nightside (right), 
showing the terms in Equations |23| and |24| Above the wind base on both sides, the escaping gas is in equilibrium, as shown by the total 
at 0. The shaded region highlights gas below the outflow base, on the nightside. 


that is ionized is low enough that recombination is too 
slow to balance ionization. Above recombination 

balances photoionization. 

4.1.2. Nightside flow: advected ions mixed with a neutral 

shadow 

As seen in Figure!^ heated dayside gas not only moves 
radially outwards, but also advects around the planet. 
The advected flow exceeds the sound speed and escape 
velocity, before leaving the box in a steady-state wind. 
The nightside flows are aligned along the anti-stellar ray, 
in part due to tidal gravity. 

Unlike the dayside flow, the nightside wind is not di¬ 
rectly driven by photoionization heating. Instead, flows 
moving around the planet converge on the nightside, 
at a stagnation point, and heat the surrounding gas to 
^ lO^K. Below the stagnation point, PdV work adds to 
the internal energy. Above the stagnation point, internal 
energy drives the outflow, which then cools by PdV ex¬ 
pansion. The interplay between PdV work and internal 
energy storage allows the gas to be in an energetic equi¬ 
librium. With the exception of localized Lya cooling near 
the stagnation point, where the gas is hot enough, other 
sources of cooling, including recombination, are negligi¬ 


ble. 

In the planet’s shadow, the gas has a sharply de¬ 
fined neutral outflow. As highlighted in Figure the 
planet’s nightside atmosphere contributes neutral gas to 
this flow, in a time-varying circulation pattern. This gas 
originates above the atmosphere’s initial outer radius of 
1.56 X 10^^cm= 1.04i?p, rather than the dayside value of 
i?p, since there is no photoionization in the shadow to 
move away material below this radius. We note that this 
is the only portion of the flow that retains memory of 
our atmospheric initial conditions. As shown in Figure 
this nightside neutral gas is also some of the hottest 
in our computational volume, due to enhanced heating 
from the converging flows. 

The circulation is mostly confined within a stagnation 
point where the advected, dayside gas converges. Pe¬ 
riodically the circulation cell, which may be unresolved 
turbulence, grows large enough so that the circulating 
neutral gas mixes with the advected ionized gas and then 
gets dragged outwards, beyond the stagnation point. 
Since the planet is optically thick to photoionization, the 
only ionized gas in this outflow is from advection. Out¬ 
side of the planet’s shadow, the gas is optically thin to 
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Neutral fraction 



Fig. 6. — Velocity vectors overplotted on the neutral fraction 
midplane slice show gas circulating on the nightside of the planet, 
with ever changing patterns of rolls, shown at 5 X 10^s (upper) and 
6 X lO^s (lower). This motion ultimately contributes neutrals to the 
outflow in the planet’s shadow. A stagnation point is visible along 
the anti-stellar ray at ~ l.Ti^p, between the inward circulating 
material and the outflowing material. 

photoiouizatiou aud, thus, it is iouized. Recombiuatiou 
iucreases the uightside ueutral fractiou, but ouly slightly, 
because the recombiuatiou timescale is loug compared to 
the time it takes gas to flow out of our box. 

4.2. Mass loss rates and low-flux results 

We calculate the iustautaueous mass loss rate by takiug 
the average of the mass flux through a spherical shell, 
with a fluite thickuess of oue cell, usiug 

= (pvr) 47rr^. (25) 

For comparisou, we also calculate the mass loss flux 
through the faces of a cube with dimeusious 2r aud ob- 
taiu comparable mass loss rates. 

As showu iu Figure we hud a steady-state mass loss 
rate of 1.9 x 10^^ g s^ for our fiducial model. It is iu- 
terestiug to uote how the mass loss rate tracks the iuput 
flux over time. Before the flux plateaus to its coustaut 
value, the mass loss rate is larger at smaller radii - iu- 



Time [10M 

Fig. 7. — Mass loss rates, calculated at the Roche lobe (2.7Rp, 
blue) and the edge of the box (5Rp, green), reach a steady-state 
value of 1.9 X 10^^ g s“^, for a model with an incident flux of 
2.02 X cm“^ s“^. The evolution of the mass loss rate tracks 
the incident UV flux (orange). 

dicative of the time it takes for the heated gas to expaud 
aud escape. The mass loss rate reaches steady-state at 
^ 3 X lO^s, just after the flux reaches its coustaut value 
of 2.02 X 10^^ cm ^ s Iu steady-state, the mass loss 
rate is the same at all radii above the Roche lobe radius 
of 2.7Rp. 

To examiue the mass loss rate’s scaliug with UV flux, 
we also carry out au ideutical simulatiou with lower iuci- 
deut flux. For a simulatiou with a maximum photou flux 
of 2.02 X 10^^ cm ^ s au order of maguitude lower 
thau our fiducial model, we hud the outflow evolves with 
a similar structu re to that of the high flux model, de¬ 
scribed iu Sectiou IrTl This lower flux simulatiou reaches 
a steady state mass loss rate of 2.2 x 10^^ g s“^. 

4.3. Agreement with ID models 

We compare a ID slice from our simulatiou to the 
model of MC09 re-calculated for our parameter values. 
The ID model of MC09 is a good poiut of comparisou 
because it also iucludes iouizatiou balauce, tidal gravity, 
heatiug aud cooliug terms. While both the day aud uight 
sides are available from our 3D model, ouly the dayside 
is available from the ID models. 

Figure shows the steady-state deusity, velocity, tem¬ 
perature, aud iouizatiou profiles from both models, aloug 
the substellar ray. We hud that the r = I surface, de- 
uoted by au x, is at —Rp iu both models. Above the 
r = I surface, the two wiud solutious are comparable. 
Both models trausitiou from subsouic to supersouic ve¬ 
locities at similar dayside locatious relative to the plauet. 
We defiue this souic poiut usiug each model’s adiabatic 
souud speed, rather thau isothermal souud speed. The 
dayside souic poiut is iuterior to the plauet’s Roche lobe 
radius (2.7Rp). Wiud temperatures from both models 
agree withiu 25%. The verificatiou of our choice of iuitial 
parameters, i.e. r{Rp) = I, aud agreemeut with compa¬ 
rable ID results based ou the models of MC09 validates 
that our model is freely aud self-cousisteutly settiug its 
parameters, rather thau haviug fixed bouudary coudi- 
tious. 

For both models, the r = I locatiou coiucides with the 
sharp trausitiou iu iouizatiou fractiou at the wiud base. 
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x/Rp 

Fig. 8. — Agreement bet ween our model along the substellar ray 
(solid) and the ID model of|Murray-Clay et al.|(|2009|) (dotted), for 
the same parameter values. Uircles indicate sonic points and x’s 
denote the r = 1 surface to photoionization. The unshaded regions 
are those that we model physically. 


In the original MC09 high flux models of a less inflated 
planet, outflows with such a rapid ionization change are 
characterized as radiative/recombination-limited. How¬ 
ever, our highly inflated planet’s outflow is energy- 
limited. The majority of energy deposited by photoion¬ 
ization is converted into PdV work. Advection, rather 
than recombination, balances photoionization at the very 
base of the wind. Our planet’s highly inflated parameters 
put its outflow near the boundary of the two regimes. 

We compare our mass loss rates to the MC09 ID 
model, rerun for our stellar and planetary parameters. 
We estimate the mass loss rate for each ID model at a 
radius r by multiplying the substellar mass loss rate by 
and a geometric correction factor, to account for 
differences in the received flux over the planet’s surface. 
For energy-limited outflows, MC09 used a constant cor¬ 
rection factor of 0.26, derived from mass loss models with 
different fluxes that showed M oc Fuy. 

Our mass loss rates for maximum photon fluxes of 
2.02 X 10^^ and 2.02 x 10^^ cm“^ s“^, respectively, agree 
with the ID model loss rates and UV flux scaling, as 
shown in FigureWe are unable to make direct compar¬ 
isons to other mass loss predictions, due to our choice of 
model parameters. However, given that MC09 had mass 
loss rates in good agreement with other ID models for 
HD 2 09458 fY^pOil |Tian et al 


2005 Garcfa Munoz 


2007), our agreement with the IViCOO model recalculated 
for our parameters suggests agreement between ours and 



Fig. 9.— Mass loss rates as a function of UV flux from this 
work (points) and the ID model of MC09 (solid line), rerun for 
the parameters used in this work. The mass loss rate scales as 
M oc -F^y for energy-limited escape. 
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Fig. 10.— Channel maps of the spatial distribution of Lyo; obscu¬ 
ration at line center and then off-center by the specified velocities. 
The upper panel shows the initial obscuration and the lower panel 
shows steady-state outflow obscuration, at 1.39 x 10®s. The dashed 
circle shows the spatial extent of the stellar disk. 

other previous ID models. 


5. PREDICTED LYMAN-ALPHA EXTINCTION 

For outflows from our pure hydrogen planet, we cal¬ 
culate the predicted obscuration of stellar Lya emission 
during transit. We focus on Lya because it has been 
observed for several planets including HD 209458b, HD 
189733b, and GJ 436b. Our calculation serves as a first 
step toward motivating future observations of a wider 
range of exoplanets, including WASP-17b. 

To determine the spatial extent of Lya absorption, we 
use the temperature and ionization fraction of each cell 
in our 3D simulation to compute a Voigt line profile, 
Doppler shifted by the cell’s b ulk radial velocity. We 
use values from Morton (2003) for the Lya oscillator 
strength, Einstein A coefficient, and line center wave¬ 
length. In this work, we calculate the obscuration pro¬ 
duced by a planet passing across the stellar disk mid¬ 
transit. We assume that the star has = Rq. 

Given that our simulation is in 3D, we can examine the 
integrated absorption spatially, in 2D. Ghannel maps of 
the obscuration, 1 — exp(—r), initially and in steady- 
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Fig. 11. — Average Lya obscuration integrated over the stellar 
disk, caused by steady-state outflows. 

state, are showu iu Figure Extiuctiou from the out¬ 
flow is axisymmetric across our simulated volume. The 
velocity depeudeuce of the extiuctiou is showu iu these 
maps aud more clearly seeu iu Figure pTj which shows the 
total obscuratiou averaged over cells iu the stellar disk. 
To make this figure cousisteut with traditioual sigu cou- 
veutious, the velocity sigus have beeu reversed aud are 
thus the opposite of other velocity values iu this paper. 
The obscuratiou is roughly symmetric about hue ceuter 
aud drops to 5% at ± ^ 50 km s“^. At higher velocities, 
iu the hue wiugs, obscuratiou of a few perceut comes ouly 
from the plauetary disk. Coufouudiug geocoroual emis- 
siou aud absorptiou by the iuterstellar medium make ob- 
servatious of Lya absorptiou challeugiug uear hue ceuter, 
where wiud obscuratiou is most promiueut. 

6. SUMMARY AND DISCUSSION 

We have simulated global hot Jupiter outflows iu 
three-dimeusious (3D). The outflow is self-cousisteutly 
lauuched by photoiouizatiou heatiug, which we simulate 
usiug a uew implemeutatiou of plauar radiative trausfer 
iu the Atheua hydrodyuamic code. Our code is compati¬ 
ble with multiple levels of static mesh refluemeut (SMR), 
distributed across parallel processors, aud is freely avail¬ 
able for further use. 

We drive outflows from a highly iuflated hydrogeu 
plauet, illumiuated by a plaue-parallel source of iouiziug 
radiatiou. We hud that photoiouizatiou heated super- 
souic outflows emerge ou the dayside aud advect arouud 
the plauet, lauuchiug uightside outflows uot captured iu 
earlier ID simulatious. Outflows are everywhere iouized, 
except iu the plauet’s uightside shadow, where the out¬ 
flow is largely ueutral. Ou the plauet’s uight side, a stag- 
uatiou poiut iu the flow separates outflowiug gas from 
a time-variable circulatiou regiou. Perhaps couuteriutu- 
itively, the louger resideuce time of gas iu this regiou al¬ 
lows it to reach temperatures higher thau those achieved 
ou the plauet’s dayside. For our fiducial parameters, gas 
withiu the uightside staguatiou regiou is the ouly portiou 
of the flow that is substautially radiatively cooled. 

For the highly iuflated plauet we cousider, we hud 
mass loss rates of 1.9 x 10^^ g s“^ aud 2.2 x lO^^g s 
for photou fluxes of 2.02 x 10^^ cm“^ s“^ aud 2.02 x 
10^^ cm“^ s“^, respectively. The outflow is margiually 
euergy-limited. Our mass loss rates are cousisteut with 


the ID escape model of MC09, reruu for the same pa¬ 
rameters. Aloug the dayside substellar ray, we also hud 
remarkable agreemeut betweeu the outflow structure iu 
our simulatious aud the ID model. 

A beuefit of our 3D model is the ability to examiue 
uot ouly day aud uight differeuces iu outflows but also 
positiou-depeudeut extiuctiou. Neutral, shadowed gas, is 
a major coutributor to the absorptiou of stellar Lya emis¬ 
sion - a key predicted observable. lutegratiug through 
the box, we hud that the Lya absorptiou of the escapiug 
gas obscures stellar emissiou out to ± ~ 50 km s“^. 

The work preseuted iu this paper is the startiug poiut 
for more realistic modeliug of atmospheric escape. While 
we have simulated asymmetries iu stellar heatiug, the 
outflow will also be shaped by the Coriolis force, mag- 
uetic fields, aud iuteractious with the stellar wiud. Iu 
additiou, true atmospheres are irradiated by a full spec¬ 
trum of euergetic photous. As our work is motivated by 
observatious of mass loss from close-iu plauets, iu the fol- 
lowiug paragraphs we review how these physical effects 
may chauge the observed trausmissiou spectrum. Siuce 
Lya observatious do uot yet exist for our plauet or the 
similar WASP-17, we discuss the impact of these physi¬ 
cal processes iu the coutext of the observed high v elocity 
Lya obscuratiou at TlQQkm s“^ for HD 209458b (Vidal- 
Madjar et al.||2003). Bear iu miud that the two planets 
have difierent parameters and thus will not have identical 
transmission spectra, even with additional physics. 

To distinguish 3D geometric effects on the outflow, we 
have used a monochromatic flux source and compared 
our results to ID models, which also used monochro¬ 
matic flux sources. Energetic effects involved in using a 
full stellar spectrum were not studied in this first 3D pho¬ 
toionization study. Due to the wavelength dependence 
of the photoiouizatiou cross section, cTph, a full stellar 
spectrum can smear out the r = 1 surface where pho¬ 
tons are absorbed, iucreasiug the thickness of the wiud 


base (jTrammell et al.||2011t |koskinen et al.||2013| |2014|) . 
These changes may increase the column density near the 
base of the wind or increase the radial extent of neu¬ 


tral hydrogen (Koskinen et al. 2013). These effects can 
contribute to a broadened line prome. At very high in¬ 
cident fluxes, photoionizing X-rays, which have a lower 
photoiouizatiou cross section per hydrogen nucleus and 
hence are deposited deeper in the atmosphere, can inject 
heat quickly enough to be a dominant driver of the out¬ 
flow ( Owen fc Jacksou||2Q12 ). The resulting wind struc¬ 
ture IS analogous to that modeled for UV-driven winds 
but begins deeper in the atmosphere, increasing the to¬ 
tal wind column. To study this effect, additional spectral 
bins including photoiouizatiou by X-rays can be included 
in future versions of our model. 

High velocity neutral gas may arise in a bow shock 
that marks the interaction between the outflow and the 


stellar wind. Tremblin & Chiang (2013) find that charge 
exchange between outflowing neutral atoms and faster 
stellar wind ions in this wind-wind interaction region can 
p roduce a sufficient fast ueutr al population. 


Trammell et al. (2011,2014) have suggested that if the 
planeUs magnetic held is dipolar, a magnetically con¬ 
fined dead zone near the equator does not participate in 
the outflow and could host an enhanced column of neu¬ 
tral gas. MC09 estimate that if the line-of-sight column 
of neutral gas were enhanced by a factor of 3-5, then 
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naturally broadened line wings from this slower-moving 
population would be sufficient to pro duce the observed 


obscuration (c.f. Figure 11 see also Ben-Jaffel 


2008). 


Regardless of wheth er magn e tic co nfinement enhances 
the neutral column, Adams (2011) demonstrates that 
for reasonable hot Jupiter magnetic field strengths, the 
magnetic field will redirect and possibly reduce the atmo¬ 
spheric outflow, affecting its 3D structure. For our fidu¬ 
cial model, we find that a local field strength B > 0.06G 
is sufficient for the magnetic pressure B‘^1 Sir to overcome 
the ram pressure of the wind at all radii. For example, 
at the sonic point of the wind located at ^2Rp, this local 
B corresponds to that from a dipo l e field with a surface 
strength of 0.5G. |Owen & Adamsj (|2014|) find using 2D 
models that increasing magnetic field strength can sup¬ 
press nightside outflows and reduce planetary mass loss 
rates. Because we have used Athena, an MHD code, for 
our simulations, it can be extended to study the effects 
of the magnetic field on planetary wind confinement. 

A line-of-sight column density enhancement could al¬ 
ternatively result from the confinement of the atmo¬ 
spheric outflow by the stellar wind. Depending on the 
level of stellar activity, pressure confinement by the stel¬ 
lar wind can channel outflowing g as into a smaller solid 
angle, increasing its local density (Stone & Proga|2009). 
In extreme cases, the flow is confined to a cylindrical re- 
gion on the planet’s night side. The Coriolis force (due to 
the planet’s orbital motion) will limit the extent of this 
high-density region by causing the outflow to curve. The 
Coriolis force bends trajectories on a timescale ^ 0“^, 
where O is the angular frequency of the planet’s orbit. 
For our outflow velocity of 'T’ ^ 20 km s“^, the length 
scale of this curvature is Tcurve ^ ^ 10^^ cm, about 

twice the Roche lobe radius. 

Finally, we note that the Coriolis force will produce 
curvature in the ne utral shadow highlighte d here and 
in the 2D model of Owen & Adams] (2014). Without 
such curvature, this large reservoir of neutral gas does 
not contribute to the line-of-sight absorption signature 
since the planet lies between it and the star. At dis¬ 
tances larger than I/curve, the neutral gas will curve out 
of the planet’s shadow, where it has the potential to sub¬ 
stantially enhance the transit absorption signature on 
the side of the planet opposite to its direction of mo¬ 
tion. Upon leaving the nightside, a neutral atom will 
ionize within ^ (T<jph)“^ ^ a few hours, sufficient time 
to travel approximately a planetary radius at 20 km s“^. 

Incorporating the Coriolis force into our simulation will 
clarify how stellar photoionization changes the outflow 
structure and allow us to produce more realistic, spa¬ 
tially resolved extinction maps and spectra. With the in¬ 
troduction of a stellar wind and magnetic fields, compar¬ 
isons can be made to 3D studies of these effects that did 
not includ e self-consistent photoionization driven wind 

launching (Cohen et al.|2Qll Trammell et al.|2011 

-*- 


sakos et ffi||2UI5[ ) . 


Mat- 


W hile the model itself can be expanded to include addi¬ 
tional physics, post-processing will also offer avenues for 
other comparisons to observation. Our hydrogen results 
can be translated into other species, which can be used to 
look for observational signatures of outflows. Given the 
3D nature of our code and the ability to track the time 
evolution of the gas, we can use a larger computational 


box to look for sightline specific effects that would ap¬ 
pear at ingress and egress. Such a model would allow us 
to make comparisons to observations with su fficient time 


resolution, including Kulow et al. (2014) and Lecavelier 
des Etangs et al. ( 2Q12| ) . With improved physics and 


ulated observational output, we will be one step closer to 
motivating and understanding observations of mass loss 
from hot Jupiters. 
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APPENDIX 

A. IONIZING RADIATIVE TRANSFER SETUP AND 
VALIDATION 


A.l. Implementation with SMR and Parallelization 

We have e xtended the ray-tracing radiative transfer 
method from Krumholz et al. (2007) to Athena grids 
with multiple levels of static mesh refinement (SMR), 
parallelized with MPI. We use a plane-parallel source of 
ionizing radiation, with rays aligned with the grid. As a 
result, we loop over cells along the direction of radiation 
propagation and iteratively carry out radiation updates. 
Radiative source terms, given on the right hand side of 
Equations!^ andare calculated and used to update the 
neutral fraction and energ y density, using the ope rator- 
split approach described in Krumholz et al. (2007). Key 


to calculating the photoionization source terms is the ion¬ 
izing flux in each cell. Eor rays passing through differ¬ 
ent levels of resolution, we communicate ionizing fluxes 
across each level and processor boundary. 

Eor a given timestep, our steps for ionizing radiative 
transfer on a grid with multiple SMR levels, from coarsest 
(level 0) to finest resolution are: 


1. Integrate the flux along rays on level 0. 

2. Compute ionization, heating, and cooling rates. 
Update the energy density and neutral fraction in 
each cell. 


3. Evaluate stopping criteria on level 0, and store the 
iteration timestep. 

4. Prolongate the flux from the current, coarse level 
(e.g. level 0) onto the next level of resolution (e.g. 
level 1), where the coarse grid first crosses into the 
finer grid. To prolongate the flux, copy the coarse 
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flux iuto the overlappiug hue cells, so that each 
coarse level ray sub-divides iuto eight rays with 
equal flux, at the uext level of resolutiou. 

5. lutegrate the flux aloug the uext level of resolutiou 
(e.g. level 1), aud repeat step 2. 

6. Coutiuue this level’s iutegratiou aud radiatiou up¬ 
date uutil the coarsest, level 0, timestep has beeu 
reached. 


7. Repeat steps 4-6 for the remaiuiug resolutiou levels. 

8. Restrict the ueutral fractiou aud euergy deusity 
from fiuer, higher resolutiou levels to coarser levels, 
e.g. level 2 to level 1, level 1 to level 0, by averag- 
iug the fiuer values of the eight overlappiug cells to 
obtaiu oue value for the coarse level ray. 


As described above, wheu SMR is used, the coarsest 
grid level is evaluated first, followed by subsequeutly fiuer 
levels. For each level, we iutegrate the flux at each cell 
usiug Equatiouj^with a discrete calculatiou of the optical 
depth, Equatiou Calculated fluxes are used to com¬ 
pute radiative source terms aud update the ueutral frac¬ 
tiou aud euergy deusity iu each cell. As these radiatiou 
updates are computatioually cheap relative to hydrody- 
uamic updates, multiple radiatiou updates cau be car¬ 
ried out before euteriug the hydrodyuamics iutegrat or. 
We use the stoppiug criteria of Krumholz et al. (2007) to 
determiue wheu the coarse grid radiatiou update should 
termiuate. To eusure that propagatiou speeds are the 
same across differeut resolutiou levels, we require fiuer 
levels to complete radiatiou updates for the same elapsed 
time as that of the coarse level. 

With SMR, the flux is ueeded at the bouudaries be- 
tweeu differeut levels of resolutiou, so we prolougate the 
radiative flux from lower resolutiou levels to cells iu over¬ 
lappiug higher level cells. Resolutiou level bouudaries 
aud overlappiug regious are determiued duriug the sim- 
ulatiou iuitializatiou. Prolougatiou occurs at every time 
step, ouly at level bouudaries, siuce the flux is theu iter¬ 
atively calculated aloug grid-aligued rays. Photou uum- 
ber is couserved across level bouudaries, siuce we copy 
the radiative flux iu each coarse cell iuto the fiuer cells, 
which have smaller areas. Recall that Atheua requires 
each level of resolutiou to be smaller by a factor of 2, iu 
each dimeusiou. 

Eollowiug the radiatiou update, we do uot restrict the 
radiative flux from the hue grids back to the coarse grids. 
Iu our mass loss simulatiou, the ambieut gas is optically 
thiu aud the plauet is optically thick, so auy flux that 
ueeds atteuuatiug will be visible ou the coarse levels as 
ueeded. We do, however, iuclude a restrictiou step from 
the fiuest levels back to the coarse levels, to recombiue 
fiuer rays together, for the hydrodyuamic variables, E 
aud pn- 

Wheu usiug parallelizatiou, extra commuuicatiou is re¬ 
quired to coordiuate, seud, aud receive iouiziug fluxes 
ou differeut processors. SMR prolougatiou aud restric¬ 
tiou operatious are complicated by the fact that values 
are distributed ou differeut processors. To appropriately 
direct fluxes to the grid levels spread across multiple 
processors, we store aud commuuicate processor ideuti- 
fiers for the uuderlyiug coarse aud overlappiug hue grids. 


Density g/cc] 



X-axis [10''^ cm] 


Fig. 12.— 2D midplane slice showing the density of an ionization 
front in a uniformly dense gas, simulated with four levels of resolu¬ 
tion. To highlight resolution level boundaries, densities in this slice 
have been rendered with different levels of transparency for each 
resolution level - coarse levels are more transparent. Radiation is 
incident from the -x direction. 

With this iuformatiou, flux prolougatiou uow occurs as 
two separate steps - a seudiug of the flux values from the 
coarse grid aud a receipt of the values ou the hue grid. 

Although parallelizatiou speeds up our simulatiou, 
there are restrictious. Eor our plaue parallel atteuuatiou 
of the flux, grids further dowustream iu the directiou of 
radiatiou propagatiou must wait for iutegratiou to hap- 
peu upstream. While radiatiou updates are computed ou 
the grids iu order aloug the directiou of radiatiou prop¬ 
agatiou, there are uo barriers to carryiug out radiatiou 
updates ou the grids distributed orthogoually to the radi¬ 
atiou propagatiou. As previously uoted, we also require 
coarse levels to complete their radiatiou updates before 
fiuer, overlappiug regious cau. 


A.2. Implementation Validation 

We examiue the accuracy of our radiatiou module with 
SMR aud parallelizatiou, by compariug our code to aua- 
lytic expectatious for a plauar iouiziug source iucideut ou 
a uuiform, optically thick ueutral medium. The evolu- 
tiou of the iouized regiou is akiu to that of au HII regiou 
Stromgreu sphere, but with a plauar geometry. 

We set up a (50 pc)^ box of ueutral gas iouized by a 
plaua r source with a photou fl u x of 3 .67 x 10^ cm“ 


.-2 g -1 


as m Geudelev & Krumholz (2012). To test our im- 
plementatiou, we use a grid with four levels of resolu¬ 
tiou, distributed amoug 8 processors, with each processor 
hostiug four grids - oue for each level of resolutiou. All 
resolutiou levels were ceutered withiu the box, with the 
coarsest three levels each haviug 32 cells, aud the fiuest 
resolutiou level haviug 48 cells. To appropriately make 
comparisous with au HII regiou, we iuitialize the gas to 
au ISM appropriate chemistry. Cousideriug atomic hy- 
drogeu with helium, we use the meau gas mass per hy- 
drogeu uucleus of /in = 2.3 x g aud the meau 

mass, /i = 2.1 X 1Q ~^^ g. We use a carb ou abuudauce 
of ac = 3 X 10“^ (|Sofia & Meyer 2001). We iuitial¬ 
ize the gas to a uuiform ueutral medium with uumber 
deusity no = 63 cm“^ aud au isothermal souud speed 
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Cs,o = 5.4 X 10^ cm s“^. Unlike our hydrogen mass loss 
simulation, we use = nH+ + pac/^ld/in), where ac 
is the carbon abundance of the gas. Since neutral car¬ 
bon has an ionization threshold below 13.6 eV, is the 
dominant ionization state in the neutral ISM, and this 
provides a minimum abundance of free electrons. 

For HII regions, different radiative source terms are 
needed i n th e flu id eq uations. We include all terms from 
Sections |2.3| and [2^ except for hya line cooling. Addi¬ 
tionally, we add contributions from collisional ionization 
to our ionization balance, usi ng the rate per unit volume 
of Tenorio-Tagle et al. (1986) 


^coll 


(5.84 X 10 ^^cm^ s n^pn 


(-13.6 eV)/kT 

K 


(Al) 

We also include the optically thin cooling terms for 
molecular cooling, line cooling f rom O, N, and Ne, and 
free-free cooling, as described in Krumholz et al. (2007). 
A full inclusion of heating and cooling terms is necessary 
to thermostat the ionized gas to 10^ K. 

We can qualitatively assess the success of our radiative 
transfer implementation by looking at a 2D slice of the 
gas density, shown in Figure Note that to make the 
different SMR regions visible, we rendered the product 
of the density and a constant transparency value, with 
greater transparency values for coarse resolutions. Radi¬ 
ation incident from the left edge of the box photoionizes 
the gas, increasing its temperature and lowering its den¬ 
sity. The ionized region is separated from the neutral 
region by an ionization front (red), which we resolve to 
one grid cell. As a result, the front appears thinner in the 
higher resolution regions, and it is also denser - to con¬ 
serve mass. The ionization structure is more clearly seen 
in Figure which gives a ID view of the density, tem¬ 
perature, neutral density, and ionization fraction. The 
ionization front is visible as the peak in density, sepa¬ 
rating the lower density ionized region from the higher 
density neutral background. 

Quantitatively, we can compare our simulation output 
to analytic calculations of ionization front expansion in 
HII regions. When an ionization source first turns on, 
the initial rate of photoionization overwhelms the rate 
of recombination. The result is a rapid phase of ioniza¬ 
tion and expansion, faster than the local sound speed. 
As this R-type expansion slows down due to reaching 
near ionization equilibrium at the Stromgren length, D- 
type expansion occurs, driven by an overpressure in the 
photoionized region. At times much greater than the 
sound crossing time of the ionized region, the photoion¬ 
ized gas will reach uniform pressure and density and exert 
a pressure on the swept up front of neutral material that 
bounds the ionized region. We can use the requirement 
from momentum conservation, that the photoionized re¬ 
gion’s force on the neutral gas be equal to the neutral 
gas’s rate of momentum change per unit area, to derive 
the size of the ionized region. The ionized region’s pres¬ 
sure on the swept up, neutral gas is 


-Pi = Cs.iPi = V (A2) 




Fig. 13.— Ionization front structure at t = 3.8 Myr, with total 
density (purple), neutral density (blue), and temperature (orange), 
shown in the upper plot. The lower plot shows the ionization frac¬ 
tion. The das hed line in the upper plot is the front location, used 
for Figure^] To the right of the front location, the neutral density 
is equal to the total density and thus hidden. 


is the mean gas mass of the ions. We assume that the 
electron number density equals the ion number den¬ 
sity rii because the photoionized region is in ionization 
equilibrium. At late times, the gas will be in thermal 
equilibrium, so we assume that is a constant. The 
time rate of change of the neutral material’s momentum 
is 

^ ^ (A3) 

where A is the cross-sectional area of the front and we 
have assumed that the distance over which mass has been 
swept. A/, is comparable t o th e siz e I o f the ionized re¬ 
gion. Equating Equations |A2| and |A3| and assuming a 
similarity solution of the form I (x yields 


l{t) = c 


yi2 V riop 


2/5 


(A4) 


where I is the size of the photoionized region, Cs^i is the Eignrep^ shows the agreement between our simulation 

sound speed in the photoionized region and pi = /iH/2 and analytic approximation for a coarse (black) and flner 
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Fig. 14. — Comparison of the ionization front location for a coarse 
domain (black), with an SMR region (blue), to D-type analytic 
approximation (red) as a function of time (top), with residuals 
(bottom). 


SMR regiou (blue). For this plot, the frout locatiou, or 
size of the iouized regiou correspouds to the x-locatiou 
of maximum deusity, for a slice through the box at fixed 
y aud z; au example of this locatiou is highlighted by the 
dashed hue iu Figure Choice of other slices orthog- 
oual to the radiatiou directiou yields similar results, as 
the gas is symmetric iu both of these directious. We de- 
termiue the correspoudiug aualytic frout locatiou, usiug 
the simulated iouized souud speed as iuput iuto Equatiou 


A4 The result we hud is that the frout locatiou agrees 


with the aualytic expressiou to withiu 5%. The oscilla¬ 
tory behavior of the residuals is a result of output that 
is uot commeusurate with the propagatiou time across 
oue cell. Nevertheless, we hud that the coarse aud hue 
regious have cousisteut agreemeut with the aualytic ap- 
proximatiou. While we have tested the physical validity 
of our radiative trausfer algorithm with SMR aud MPI, 
we have uot examiued its efficieucy aud scaliug. 


B. CONVERGENCE OF MASS LOSS SIMULATIONS 

To validate our uumerical results, we ruu a couver- 
geuce test by decreasiug the resolutiou of the plauet’s 


atmosphere by a factor of two. Siuce additioual refiue- 
meut levels iucrease the resolutiou of our plauet, uot the 
eutire box, we re-ruu our fiducial simulatiou with four 
levels of refiuemeut, iustead of five. 

As highlighted iu the ID profiles iu Figure these 
two simulatious are couverged aud produce comparable 
outflows. As expected, differeuces arise ouly iu the sta- 



x/Rn 



x/Rn 



X/Rn 


Fig. 15. — High resolution (solid) and lower resolution (dashed 
blue) simulations have well-converged outflow results, as shown in 
these ID profiles. 


bility of the static atmosphere, siuce lower resolutiou de¬ 
creased the uumber of resolved scale leugths. Thus, the 
ouly differeuces iu this resolutiou test were seeu arouud 
the edge of the atmosphere, uot iu the escapiug gas. We 
hud uo differeuce iu the calculated mass loss rates for 
the two resolutious, further suggestiug that our results 
are couverged. 

Iu the maiu text, we preseut our high resolutiou results 
for our fiducial parameters. Due to the couvergeuce we 
describe here, from our resolutiou tests of four aud five 
levels of refiuemeut, we use ouly four levels of refiuemeut 
for the low-fiux simulatiou described iu the maiu text. 
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